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Phase-space representations based on coherent states (P, Q, Wigner) have been successful in the 
creation of stochastic differential equations (SDEs) for the efficient stochastic simulation of high 
dimensional quantum systems. However many problems using these techniques remain intractable 
over long integrations times. We present a number-phase Wigner representation that can be unrav- 
eled into SDEs. We demonstrate convergence to the correct solution for an anharmonic oscillator 
with small dampening for significantly longer than other phase space representations. This process 
requires an effective sampling of a non-classical probability distribution. We describe and demon- 
strate a method of achieving this sampling using stochastic weights. 

PACS numbers: 42.50.Lc,02.70.-c,02.50.Ey,02.50.Cw,03.75.Kk 



I. INTRODUCTION 

Coherent states have traditionally been used as the ba- 
sis for the creation of phase-space representations [2]. A 
major use of these representations has been the produc- 
tion of stochastic differential equations that can greatly 
reduce the dimensionality of a problem. For example, a 
set of iV harmonic oscillators which would typically re- 
quire a density matrix with D"^^ components (where D 
is the number of elements you have in your truncated 
basis) to solve directly can be changed to a set of only N 
stochastic differential equations. Not only are the mem- 
ory requirements of the generated stochastic differential 
equations logarithmically smaller, they also typically in- 
tegrate much faster than a direct master equation ap- 
proach [nnigiiigiH]. 

Unfortunately, some problems have severely limited in- 
tegration times or cannot be solved at all using standard 
phase-space methods [TJ [13 . A coherent state ba- 
sis is inappropriate for many problems and we demon- 
strate that there is a class of quantum systems for which 
a number-phase based alternative will do better. A va- 
riety of number-phase 'Wigner-like' representations have 
been created in the past. The primary reason for such a 
diversity of solutions to a problem, which is quite clear 
in the coherent state case, is the ambiguity of the phase 
operator in quantum mechanics. The problem was origi- 
nally investigated by Dirac when he attempted to canon- 
ically quantise the electric field using number and phase 
as canonical operators, however it was soon found that 
such an approach leads to irreparable contradictions |llj . 
A more algebraic approach was taken by Susskind and 
Glowger, which resulted in a phase operator that had 
many of the required properties but was not Hcrmitian 
[12j . This issue was investigated and improved using a va- 
riety of techniques created by Pegg and Barnett [HI [14], 
and many number-phase-space representations were cre- 



ated using their phase operators [151 HI] • More recently 
Moya-Cessa has presented an number-phase Wigner rep- 
resentation based on the Susskind and Glogower operator 
directly [T7] which does not require an extention or trun- 
cation of the fock space like Pegg and Barnett 's phase 
operators. However, all these approaches were primar- 
ily motivated with visualising quantum states and lacked 
the operator correspondences required for the generation 
of stochastic differential equations. This paper presents 
a number-phase Wigner representation that can be used 
to generate stochastic differential equations. We inves- 
tigate the properties of this distribution and use it to 
demonstrate efficient numerical results. 

We find that a coherent state in the number-phase 
Wigner representation has a non-classical probability dis- 
tribution, where sections of the distribution are negative. 
Simulation of such states is important for comparative 
and practical applications of the representation. Thus, 
we develop a method for sampling non-classical distribu- 
tions using stochastic weights. Stochastic weights have 
been used previously to improve convergence of positive 
P simulations [18 and in the simulation of stochastic 
conditional master equations [19] . 



II. NUMBER-PHASE WIGNER 
REPRESENTATION 

We define our representation in such a way that a mas- 
ter equation can be converted directly to the equation of 
motion for the phase-space function, similar to Wigner 's 
original formulation [20] , as opposed to defining it by the 
ordering of the operators. We take the definition: 
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where n G 0,1/2,1,3/2,2..., (j) G (0,27r), p is a density 
matrix and \p) is a Fock state {p € 0,1,2...) the eigen- 
state of the number operator — p\p). The sum in- 
creases by increments of 1 (not 1/2), ensuring that the 
Fock states are always integer numbers. While in many 
previous phase-space representations the domain of the 
variables is the set of eigenvalues of particular operators, 
in this case we require n to take half integer values as a 
convenient way to ensure that the representation is com- 
plete. Consider 
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thus n and to must take on half integer values to ensure 
that that all density matrix elements (such as (n|/5|n-|-l)) 
can be extracted. We can see that this definition has 
many of the properties expected from a Wigner repre- 
sentation: 

i) W{n, (j)) is always real 
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ii) Integrating out the phase variable return the expected 
distribution for number 
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d(j) VF(n, 0) = (n\p\n) = P{n). 



(4) 



iii) Summing out the number variable returns the ex- 
pected distribution for phase 
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where a,b € 0,1,2... and {(p) is the eigenfunction of 
gi$ _ Y^'^^Q \n){n + 1| the exponential phase operator 
as defined by Susskind and Glogower [T^ . 

The main advantage of this representation is that it is 
straightforward to derive the operator correspondences 
between a master equation and a partial differential equa- 
tion for the evolution of the number-phase representa- 
tion. For example consider the correspondence for the 



number operator, replacing pn into our definition ([T]) 
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Using similar techniques for the other operators we find: 
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The expectation values in the phase operator correspon- 
dences can be re-written as correlations with the vacuum 
as /q #'e-2™'^V(n, 0') = (0|/5|2n). Thus for large 
fields where occupation of the vacuum is low, we can 
neglect these terms and our operator correspondences be- 
come analogous to the coherent Wigner case where x n 
and p — > (a; and p are the position and momentum op- 
erators receptively). The canonical commutator and the 
Hermiticity of the phase operator is restored in the large 
field limit so this result is not surprising. 

As one must sample the distribution before any numer- 
ical simulation, it is useful to examine the distribution for 
some typical physical states. First let's consider a num- 
ber state p = |TO)(m| 
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The state is completely defined in number but unknown 
in phase as one would expect. Next let's consider a co- 
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Finally, we also need to determine the expectation value 
correspondence. Investigating number and phase oper- 
ator ordering is probably of little practical use as rele- 
vant observables are typically not a simple function of 
ordered number and phase operators, and re-arranging 
such functions into a specific ordering is not simple. For- 
tunately, we can find a connection between moments of 
our Wigner function and anti-normally ordered ladder 
operators. Consider 
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are kept strictly positive and are a measure of the rel- 
ative confidence of the data [211 l22]- However, in re- 
cent times weights have become more common as a use- 
ful tool in a variety of fields: in the simulation of con- 
ditional master equations [19 , and to improve conver- 
gence in stochastic simulations of multimode bosonic and 
fermionic fields [531 HI]- In some recent applications of 
weights it has been advantageous to allow weights to take 
negative values, whether this is because of the weights 
become negative during the implementation of an algo- 
rithm like in the field of engineering [251 HSl 127] or for 
sampling purposes in monte carlo simulation of fermionic 
systems [28t[^. 

We present two methods based on weights. The first 
is a probabilistic method that we will prove is the op- 
timal way of sampling a non-classical probability distri- 
bution in general. The second method is a deterministic 
method, which may have an advantage when the compu- 
tational requirements of the first method are too high for 
particular distributions. We also quantify how these two 
methods can be compared. 

We define weighted averages as follows 
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Note we made a substitution on the 4th line and let n 
n + ^ xhe identity above is the most straightforward 
way to make expectation value correspondences for most 
physically interesting observables. 



III. NON-CLASSICAL PROBABILITY 
SAMPLING 

In section |lT] we have found that a coherent state in 
the number-phase Wigner representation does have some 
negativity. Only strictly positive probability distribu- 
tions can be sampled using traditional sampling tech- 
niques, so it is necessary to introduce an additional de- 
gree of freedom to allow unorthodox averaging. We do 
this by introducing a weight for each path, and in par- 
ticular we allow these weights to be negative. Tradi- 
tionally weights have been used when combining data 
sets of varying confidence, for example combining the re- 
sults of different experiments where a different number 
of measurements have been taken. In this case, weights 
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We begin with the probabilistic method. Consider a 
probability P(x) that is normalised, but not strictly pos- 
itive. We assume -P(x) can be split into a strictly positive 
function, P~''(x) > V x, and strictly negative function, 
P~(x) < V X, as follows 



P(x) = P+(x) + P"(x). 



(12) 



We further require that the functions -P*(x) have fi- 
nite integral, specifically — J dx P^(x) (note that 
C"*" -I- C~ = 1). We can then define two strictly pos- 
itive normalised probability distributions as -P*(x) = 
(x) /C^ , which can be sampled in a traditional man- 
ner. Now let us consider the average produced if we 
sample the random variables xf^ from the distributions 
P*(x) with rates A*, and assign each path a weight 
u!± — I A* for a total of N samples where 1 > A* > 
and A+ A^ = 1. For a sufficiently large N we can as- 
sume A^A+ paths will be sampled from P~''(x) and A^A^ 
paths will be sampled from P~ (x) , the average produced 
by this technique, using ( |lip , is 

Er"-+/(x+)+Er^^^/(xr) 
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In the limit of large N the sum's of random variables in 
( 13 1 converge to the integral of the probability distribu- 



tion, as required by the definition of a random variable. 
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thus we can replace 

7(x) = /dx/(x)P+(x)C++ / dx/(x)p-(x)C- 



I rfx/(x)(P+(x)+p-(x)) 



/(x) 



dx/(x)P(x). 
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Thus, the non-classical distribution P(x) is sampled 
by using the dual distributions and probabilities out- 
lined above, and computing averages using the weighted 
method. The weights in this method are normalised in 
the sense that E[u!] = 1, but this can be scaled arbitrarily 
if desired. There is still a large degree of freedom in what 
we choose for and , which needs to be further in- 
vestigated to determine the optimal sampling technique. 
Optimisation requires a measure of the accuracy of the 
calculated means, which is traditionally measured using 
the variance. This approach is validated by the central 
limit theorem, which can be extended to apply to the 
weighted mean case. Assuming the weights themselves 
can be treated as random variables, and remembering to 
take into account correlations, we find: 
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The variance reduces in the limit of large sampling, but 
to make the most efficient use of resources, we want to 
minimise the variance for a fixed N. As mentioned above, 
we have a degree of freedom in our choice of both and 
which can be optimised. The terms that depend on 
/(x) will have a convoluted dependence on P(x) and on 
the particular observable that is being calculated. The 
best we can hope achieve is to minimise E[a;2]/E[cLi]2, this 
will reduce the uncertainty in observables as a whole. Al- 
though we acknowledge that for specific observables and 
P(x)'s this is not necessarily true, from now on we re- 
strict our analysis to minimising E[w2]/E[u;]2 with respect 
to A+ and C'^. Taking the derivative of 



(c+)2 (c-)2 _ (c+)2 {i-c+y 



This shows that the optimal sampling technique is 
achieved by splitting P(x) = P^(x) -I- P~(x) such that 
(C+ — C^)2 is minimised. This can be achieved by defin- 
ing P^{x) piecewise, bounded by the points where P(x) 
crosses the x axes. Then the random variables xf^ are 
sampled from the distribution P*(x) with a rate A* = 
and assign each path a weight llj± — (C+ — C^) 



C+-C- 
for a total of N samples. 

Of course, if the form of the resulting functions is 
computationally difficult to sample, a method with sub- 
optimal sampling might be more efficient overall. An 
example of this is a simulation where we wish to sam- 
ple the coherent state from equation Splitting the 
function into two as outlined above and analytically in- 
verting the parts for sampling is not possible. While 
numerically inverting the equations would be an option, 
this would restrict our random samples to starting on a 
coarse grid, and the majority of the advantages outlined 
above would be lost. Instead we can use a simple deter- 
ministic technique to approximate the initial sampling. 
The additional speed in using this technique allows us 
to use a finer grid, which will make this technique easier 
and competitive. 

Consider sampling the function P(x), assume that the 
function is zero (or sufficiently close to zero) outside of a 
■volume unit of V. Divide this volume into equally spaced 
grid of N points, then at each point assign the random 
variable x^ the value of the point at that grid and give 
that path a weight of uii = VP{'Xi). Using (11) we see 



the average calculated using this technique gives 
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where Ax = V/N. In the limit of large N we can assume 
that the sums approximate the integrals outlined below 
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as required. We can investigate the contribution to the 
variance of a calculated mean 



EM2 A+ A- A+ 1 - A+ ' ' ' 

we find that E[lj^] /E[u;]^ is minimised when A+ = C+ or = Y. ^Ax(P(x,))^ J rfx(P(x))2 (20) 

A+ = ^Jrl , = 7^#^. The first case A+ = C+ is only 



2C+-1 — c+-c~ ■ "^^^ ^'^^^ '^^^ ~ only 
possible when C+ = 1 and C~ = as both < A"*" < 1, 
this scenario only occurs when the probability distribu- 
tion P{x) has no negativity and the weighted sampling 
technique reduces to its traditional counterpart. The sec- 
ond case is of more practical use, remembering > 1 
and < the expression for A+ = 

A 

can be sampled traditionally. Replacing these into ( 16 ) 
we obtain 



imphes 

q+'Iq- which both satisfy < A* < 1 thus they 
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Unlike the previous technique we cannot optimise this ex- 
pression, as V is fixed by the distribution itself. Typically 
the expression above will be significantly larger than the 
maximum efficiency achieved using the first technique. 
When considering a problem one needs to compare the 
ratio E[a;2]/E[ci;] for each of the two methods above in 
order to decide whether the computational difficulty im- 
plementing the first method is worth the improvement 
to the variance. The determinisitic method, if computed 
as an initial condition before simulating, can also be im- 
proved by using the breeding algorithm outlined in paper 
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IV. NUMERICAL SIMULATIONS 

We now demonstrate the effectiveness of stochastic 
methods based on our number-phase Wigner representa- 
tion by simulating a damped anharmonic oscillator. An 
analytic solution has been calculated for this problem 
with a coherent state as an initial condition in [30], and 
will be used to compare the accuracy of methods. We 
will show that in the case of weak damping, the number- 
phase Wigner representation converges over a signifi- 
cantly longer time interval than both gauge positive-P 
and truncated Wigner methods. We chose to model an 
anharmonic oscillator because it is algrebraically anal- 
ogous to simulating the nonlinear term present in the 
hamiltonian of a Bose Einstien Condensate. This nonlin- 
ear term is typically what limits the convergence of previ- 
ously developed scalable methods. The addition of damp- 
ening was neccesary because the nonlinear term is ana- 
lytic for the number-phase Wigner representation, thus in 
the undampened case it would trivially win against other 
scalable methods. To make the comparison fair we added 
dampening which is both physically reasonable and non- 
trivial for the number-phase Wigner representation. 

Consider the master equation of an anhamonic oscilla- 
tor with damping 



For a sufficiently smooth function, we can assume the 



(21) 



Applying the correspondences Q we find the number- 
phase space evolution to be 

dtW{n,(l)) — dtf,xnW{n,(j)) ~ jnW{n,(j)) 



+7y (n + 1)2 + -dl W{n + 1, 0).(22) 

Unraveling this equation into a set of stochastic differen- 
tial equations requires the application of two approxima- 
tions. First, we expand the square root 



(-+1)^ + ^5| = (" + 1) + ^ 



O(a^), (23) 



and we ignore terms of order 0{d^) and higher. This 
approximation improves in the limit of large fields. Under 
this approximation, equation ( 22 1 becomes 



8(n-f 1) 

7((n + l)W{n +1,(1))- nW{n, (/))). (24) 
We are unable to stochastically unravel this equation due 



to the second term in (24), a nonlocal differential, so we 



make the following addition and subtraction to ( 24 1 
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7((n 4- l)W{n + 1,(1))- nW{n, (j))) 
^^^^{W{n + 1,(1))- W{n, ,/))).(25) 



term on the last line of (25 1 will be small. For the coher- 
ent state, which we intend to use as an initial condition 
for our problem, this is quite a reasonable approximation. 
Ignoring this term, we make a second approximation re- 
ducing (251 to 



dtW{n, (j)) « d^xn-W^n, (j)) + 



-W{n,(t)) 



n+l) 

+-l{{n + l)W{n +1,0)- nW{n, (j))\2Q) 

We can now unravel this equation into a set of stochastic 
stochastic differential equations. With some minor in- 
spection one can recognise that the first term and second 
terms correspond to drift and diffusion in (f) respectively, 
and the last two terms correspond to a Poissonian jump 
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(27) 



Here dWt is an Ito Wiener process and dNt is a Pois- 
son jump process with a rate E[diVj] = 7E[n]di. We 
integrate ( 27 1 then compare the expectation value of po- 



sition to two other leading scaleable stochastic differen- 
tial methods: Gauge-P"'" [T5J |3T] and truncated Wigner 
[6J. Gauge-P+ is the longest-converging exact stochas- 
tic method for an anharmonic oscillator, for which it can 
generate the following equations of motion: 



dil)t 
d0t 



(x(n(l -t)- 2Gi(t + t)) - ^(1 + z)) dt + y^dW}, 
(x{n*{l -i)- 2G2(T - i)) - ^(1 + i)) dt + ^dWl 
-2T{G\ + Gl)dt + V2{G2dW^ - GidW^). (28) 



where a — e'^T^-''^* and (3 — e^^^^"^^ are the stochas- 
tic coherent amplitude and complex conjugate. These 
are complex conjugate at the beginning but may di- 
verge afterwards. We also have n = a(3* , Gi — 
0.005[i(n* -n)/2- {n* +n)/2+ \a\^] , G2 = 0.005[i(n* - 
f], and T 



n)/2 + {n + n*)/2 - ], and T = tan(6't). Position 
is calculated using x — -^E[^[uj{(3* + a)]]/E[uj] where 

m = 2e^^"^ lcos(^t). Other correspondences and more 
details are presented in [TB]. The next stochastic dif- 
ferential equation technique we consider is the coherent 
truncated Wigner method. It is an approximate method 
that assumes higher order derivates beyond diffusion can 
be safely ignored. It has a typically longer numerical 
convergence time than Gauge however, due to the 
uncontrolled approximation, there is no indication when 
the numerical solution begins to diverge from the true 
solution. The stochastic differential equations of motion 
are: 



da — I ixoi 
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FIG. 1: (Color online) Position vs. time for a dampened 
anharmonic oscillator (x = 1 a-nd 7 = 0.001) integrated us- 
ing the number-phase Wigner representation in red (dashed), 
truncated Wigner representation in green (dot dashed), gauge 
P+ in blue (dot dot dashed) and for comparison the analytic 
solution in black (solid). Sampling uncertainties are plotted 
with dotted lines in the respective colour of the numerical 
solution. Inset (a) shows the short time period gauge 
matches the analytic solution until divergent infinite paths 
dominate the evolution. The coherent Wigner representa- 
tion converges for the first section of de-phasing in inset (a), 
however inset (b) shows that it does not converge during the 
dampened ressurection of the coherent state. The method 
based on the number-phase Wigner representation converges 
over the entire integration period. The numerical integra- 
tion was performed by using the open source software package 
xpdeint, which is a new version of the xmds package [32] 

The observable is given by a; = ^E[a -I- a*]. 

The simulations using the different integration meth- 
ods are presented in Fig. [l] As we can see, the 
number-phase Wigner representation converges for sig- 
nificantly longer than the other stochastic differential 
equations and reconstructs the dampened revival of the 
coherent states successfully. This suggests the number- 



phase Wigner representation has the possibility to inte- 
grate previously difficult quantum problems over signifi- 
cantly longer integration times than other coherent based 
stochastic differential equations. The anharmonic oscil- 
lator example is of particular importance to the field of 
Bose-Einstein condensation as the anharmonic term has 
analogous algebraic properties to the dominant nonlinear 
term present in the many-body master equation. 



V. CONCLUSION 

We presented a number-phase Wigner representation 
that is suitable for the investigation of quantum dy- 
namics in phase-space. Direct operator correspondences 
were derived, allowing a master equation to be written 
as partial differential equations for the number-phase 
quasi-probability distribution. Under some approxima- 
tions, these equations can be unravelled into a set of 
low-dimensional stochastic equations, hence reducing the 
dimensionality of the system. This stochastic method 
was used to model the damped anharmonic oscillator, 
and it was shown that this method converged dramati- 
cally longer than the truncated Wigner and Gauge-P"*" 
methods. This result therefore suggests that a stochas- 
tic method based on our number-phase Wigner repre- 
sentation may be of use for currently difficult problems 
where number-conserving nonlinear terms dominate the 
dynamics. 



VI. ACKNOWLEDGEMENTS 

This research was supported under Australian Re- 
search Council's Discovery Projects funding scheme 
(project number DP0556073) and the Australian Re- 
search Council Centre of Excellence for Quantum- Atom 
Optics (ACQ AO). We acknowledge the use of CPU time 
at the National Computational Infrastructure National 
Facility and thank Graham Dennis for his help with sim- 
ulations. 



[1] C. Gardiner, Quantum Noise (Springer- Verlag, 1991). 
[2] D. Walls and G. J. Milburn, Quantum Optics (Springer- 
Verlag, 2008). 

[3] C. Gardiner, Handbook of Stochastic Methods (Springer- 
Verlag, 1983). 

[4] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, 

S. M. Tan, M. J. CoUett, D. F. Walls, and R. Graham, 

Phys. Rev. A 58, 4824 (1998). 
[5] J. J. Hope, Phys. Rev. A 64, 053608 (2001). 
[6] A. A. Norrie, R. J. Ballagh, and C. W. Gardiner, Phys. 

Rev. Lett. 94, 040401 (2005). 
[7] M. T. Johnsson and S. A. Haine, Physical Review Letters 

99, 010401 (pages 4) (2007). 



[8] R. G. Ball, L. J. Byron, A. G. Truscott, G. R. Dennis, 
M. T. Johnsson, and J. J. Hope, Physical Review A 79, 
011601 (pages 4) (2009). 
[9] A. Gilchrist, C. W. Gardiner, and P. D. Drummond, 
Phys. Rev. A 55, 3014 (1997). 

[10] A. Sinatra, C. Lobo, and Y. Castin, Journal of Physics B; 
Atomic, Molecular and Optical Physics 35, 3599 (2002). 

[11] P. A. M. Dirac, Proceedings of the Royal Society of Lon- 
don. Series A 114, 243 (1927). 

[12] L. Susskind and J. Glogower, Physics 1, 49 (1964). 

[13] S. M. Barnett and D. T. Pegg, Journal of Physics A: 
Mathematical and General 19, 3849 (1986). 

[14] D. T. Pegg and S. M. Barnett, Phys. Rev. A 39, 1665 



7 



(1989) . 

[15] J. Vaccaro, Phys. Rev. A 52, 3474 (1995). 

[16] J. A. Vaccaro and D. T. Pegg, Phys. Rev. A 41, 5156 

(1990) . 

[17] H. Moya-Cessa, Journal of Optics B: Quantum and Semi- 
classical Optics 5, S339 (2003). 

[18] P. Deuar and P. D. Drummond, Computer Physics Com- 
munications 142, 442 (2001), ISSN 0010-4655. 

[19] M. R. Hush, A. R. R. Carvalho, and J. J. Hope, Physical 
Review A 80, 013606 (pages 5) (2009). 

[20] E. Wigner, Phys. Rev. 40, 749 (1932). 

[21] W. G. Cochran, Supplement to the Journal of the Royal 
Statistical Society 4, 102 (1937), ISSN 14666162. 

[22] P. Meier, Biometrics 9, 59 (1953), ISSN 0006341X. 

[23] P. Dcuar and P. D. Drummond, Phys. Rev. A 66, 033812 
(2002). 

[24] J. F. Corney and P. D. Drummond, Phys. Rev. Lett. 93, 
260401 (2004). 



[25] J. Shi, C. Hu, and C. Shu, Journal of Computational 
Physics 175, 108 (2002). 

[26] G. Arce and J. Parcdes, Signal Processing, IEEE Trans- 
actions on 48, 768 (2000). 

[27] J. Fakcharoenphol and S. Rao, Journal of Computer and 
System Sciences 72, 868 (2006), special Issue on FOGS 
2001. 

[28] M. Nedjalkov, H. Kosina, and S. Sclbcrhcrr, Large-Scale 
Scientific Computing pp. 2048 2048 (2004). 

[29] S. Y. Chang, V. R. Pandharipandc, J. Carlson, and K. E. 
Schmidt, Phys. Rev. A 70, 043602 (2004). 

[30] G. J. Milburn and C. A. Holmes, Phys. Rev. Lett. 56, 
2237 (1986). 

[31] L. I. Phmak, M. Flcischhauer, M. K. Olsen, and M. J. 

CoUett, Phys. Rev. A 67, 013812 (2003). 
[32] Project website http://www.xmds.org. 



